﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Windows.Media.Media3D;
using System.Windows.Media;
using DataSimDLL;

namespace BDSim
{
    public class Base3DModel
    {
        /// <summary>
        /// 3D模型的初始位置
        /// </summary>
        protected Vector3D InitPosition = new Vector3D(0, 0, 0);

        /// <summary>
        /// 3D模型的当前位置
        /// </summary>
        protected Vector3D NewPosition = new Vector3D();

        /// <summary>
        /// 3D模型需要平移的向量
        /// </summary>
        protected TranslateTransform3D Translate3D = new TranslateTransform3D(0, 0, 0);

        /// <summary>
        /// 真实地球与当前地球的比例
        /// </summary>
        protected double earthRadius = 6378137 * 2;

        /// <summary>
        /// 待添加的3D模型
        /// </summary>
        public ModelUIElement3D _3DModel = new ModelUIElement3D();

        public void Translation()
        {
            //暂时删除平移量
            (_3DModel.Model.Transform as Transform3DGroup).Children.Remove(Translate3D);

            //使用当前位置减去初始位置得到平移向量
            Vector3D translate = NewPosition - InitPosition;
            Translate3D = new TranslateTransform3D(translate);

            //重新添加平移量，把3D模型移动到新的位置
            (_3DModel.Model.Transform as Transform3DGroup).Children.Add(Translate3D);
        }

        public virtual void Generate3DModel(){}

        public virtual void Generate3DTrack(){}
    }

    public class Satellite3DModel : Base3DModel
    {
        ESatelliteShape _Shape = ESatelliteShape.SimpleSatellite;

        public ESatelliteShape Shape
        {
            get { return _Shape; }
            set { _Shape = value; }
        }

        List<Brush> BrushList = new List<Brush>() {{Brushes.Red}, {Brushes.Green}, {Brushes.Blue}, {Brushes.Yellow} };

        Orbit3DParam OrbitParam = new Orbit3DParam();

        /// <summary>
        /// 卫星轨道模型
        /// </summary>
        public ModelUIElement3D _3DOrbitModel = new ModelUIElement3D();

        /// <summary>
        /// 生成卫星3D模型
        /// </summary>
        public override void Generate3DModel()
        {
            GeometryModel3D geometryModel3D = new GeometryModel3D();
            Geometry3D geometry3D = new MeshGeometry3D();
            switch (Shape)
            {
                case ESatelliteShape.Sphere:
                    SphereMeshGenerator sphere = new SphereMeshGenerator() { Radius = 0.05, Center = new Point3D(InitPosition.X, InitPosition.Y, InitPosition.Z) };
                    geometry3D = sphere.Geometry;
                    break;
                case ESatelliteShape.SimpleSatellite:
                    SimpleSatellite simpleSatellite = new SimpleSatellite() { Center = new Point3D(InitPosition.X, InitPosition.Y, InitPosition.Z) };
                    geometry3D = simpleSatellite.Geometry;
                    break;
                default:
                    break;
            }
            geometryModel3D.Geometry = geometry3D;
            geometryModel3D.Material = new DiffuseMaterial(Brushes.LightGray);
            geometryModel3D.BackMaterial = new DiffuseMaterial(Brushes.LightGray);
            Transform3DGroup transform3DGroup = new Transform3DGroup();
            transform3DGroup.Children.Add(Translate3D);
            //RotateTransform3D rotate3D1 = new RotateTransform3D();
            //rotate3D1.Rotation = new AxisAngleRotation3D() {Axis = new Vector3D(0,1,0), Angle = 30};
            //transform3DGroup.Children.Add(rotate3D1);

            geometryModel3D.Transform = transform3DGroup;
            _3DModel.Model = geometryModel3D;
        }

        /// <summary>
        /// 生成卫星轨迹（轨道）
        /// </summary>
        public override void Generate3DTrack()
        {
            //轨道半径
            GeometryModel3D geometryModel3D = new GeometryModel3D();
            Geometry3D geometry3D = new MeshGeometry3D();
            
            //建立轨道模型
            SphereRingMeshGenerator orbitMesh = new SphereRingMeshGenerator() { RingRadius = OrbitParam.Radius, Center = new Point3D(0, 0, 0) };
            geometry3D = orbitMesh.Geometry;
            geometryModel3D.Geometry = geometry3D;
            //定义表面和背面材质
            geometryModel3D.Material = new DiffuseMaterial(BrushList[(int)OrbitParam.Omega/90]);
            //geometryModel3D.BackMaterial = new DiffuseMaterial(Brushes.AliceBlue);
            //定义转换组
            Transform3DGroup transform3DGroup = new Transform3DGroup();

            RotateTransform3D zRotate3D = new RotateTransform3D();
            zRotate3D.Rotation = new AxisAngleRotation3D() { Axis = new Vector3D(0, 0, 1), Angle = OrbitParam.Incline };//设置轨道倾角
            transform3DGroup.Children.Add(zRotate3D);

            RotateTransform3D yRotate3D1 = new RotateTransform3D();
            yRotate3D1.Rotation = new AxisAngleRotation3D() { Axis = new Vector3D(0, 1, 0), Angle = OrbitParam.Omega };//绕y轴旋转设置升交点经度
            transform3DGroup.Children.Add(yRotate3D1);

            
            geometryModel3D.Transform = transform3DGroup;
            _3DOrbitModel.Model = geometryModel3D;
        }

        /// <summary>
        /// 此处使用升交点赤经解算经度（注：用UTC代替了UT1，使用格林尼治平恒星时代替了真恒星时，需要再修改）
        /// </summary>
        /// <param name="TrueAnomaly">真近点角，单位为度</param>
        /// <param name="raan">升交点赤经，单位为度</param>
        /// <param name="a">半长轴，单位为m</param>
        /// <param name="UT1">儒略日，JD</param>
        /// <returns>度</returns>
        public double RAAN2Lon(double TrueAnomaly, double raan, double a, double UT1)
        {
            //TrueAnomaly = TrueAnomaly * Math.PI / 180;//度转弧度
            //raan = raan * Math.PI / 180;//度转弧度

            double GM = 3.986004418E14;
            //double We = 7.299115e-5;
            double dt = TrueAnomaly * Math.Sqrt(Math.Pow(a, 3) / GM);
            double du = UT1 - 2451545 - dt / 86400.0;
            double Tu = du / 36525;
            double Tu2 = Tu * Tu;
            double Tu3 = Tu * Tu2;
            double GMST = TimeAng(0, 0, 67310.54841) + TimeAng(876600, 0, 8640184.812866) * Tu +
                    TimeAng(0, 0, 0.093104) * Tu2 - TimeAng(0, 0, 6.2E-6) * Tu3;

            //double RAAN = Lon + GMST;
            //RAAN = RAAN % (2 * Math.PI);
            //RAAN = RAAN * 180 / Math.PI;//弧度转度
            //return RAAN;
            
            double lon = raan - GMST;
            lon = lon % (2 * Math.PI);
            lon = lon * 180 / Math.PI;//弧度转度
            return lon;
        }

        public double TimeAng(double h, double m, double s)
        {
            double result;
            result = (((h) + (m) / 60.0 + (s) / 3600.0) * 15.0 * Math.PI / 180);
            return result;
        }

        /// <summary>
        /// 设置3D模型的初始状态，参数为卫星的ECF位置及时间
        /// </summary>
        /// <param name="satECFPV"></param>
        public void SetInitPos(NodePosVel satECI)
        {
            Orbit6Element sixElement;
            TimeSpaceDll.ECI2Element(satECI, out sixElement);
            InitPosition = new Vector3D(satECI.Y / earthRadius, satECI.Z / earthRadius, satECI.X / earthRadius);
            OrbitParam.Incline = sixElement.i;
            OrbitParam.Radius = Math.Sqrt(satECI.X * satECI.X + satECI.Y * satECI.Y + satECI.Z * satECI.Z) / earthRadius;
            
            //直接用升交点确定轨道
            OrbitParam.Omega = sixElement.O;
            MainDispatcher.MainWindowDispatcher.BeginInvoke(new Action(() =>
            {
                Generate3DModel();
                Generate3DTrack();
            }));
        }

        /// <summary>
        /// 更新3D模型的的位置
        /// </summary>
        /// <param name="satECI"></param>
        public void UpDate3DModelPos(NodePosVel satECI)
        {
            NewPosition = new Vector3D(satECI.Y / earthRadius, satECI.Z / earthRadius, satECI.X / earthRadius);
            Translation();
        }

        /// <summary>
        /// 求平面法向量
        /// </summary>
        /// <returns></returns>
        double[] FaXiangLiang(double Px1, double Py1, double Pz1,
            double Px2, double Py2, double Pz2, double Px3, double Py3, double Pz3)
        {
            double[] newJiaoPoint = new double[3];
            //平面方程Ax+BY+CZ+d=0 行列式计算
            double A = Py1 * Pz2 + Py2 * Pz3 + Py3 * Pz1 - Py1 * Pz3 - Py2 * Pz1 - Py3 * Pz2;
            double B = -(Px1 * Pz2 + Px2 * Pz3 + Px3 * Pz1 - Px3 * Pz2 - Px2 * Pz1 - Px1 * Pz3);
            double C = Px1 * Py2 + Px2 * Py3 + Px3 * Py1 - Px1 * Py3 - Px2 * Py1 - Px3 * Py2;
            double D = -(Px1 * Py2 * Pz3 + Px2 * Py3 * Pz1 + Px3 * Py1 * Pz2 - Px1 * Py3 * Pz2 - Px2 * Py1 * Pz3 - Px3 * Py2 * Pz1);
            double E = Math.Sqrt(A * A + B * B + C * C);
            newJiaoPoint[0] = A / E;
            newJiaoPoint[1] = B / E;
            newJiaoPoint[2] = C / E;
            return (newJiaoPoint);
        }
        /// <summary>
        /// 求平面法向量
        /// </summary>
        public double GetLongitude(Point3D p1, Point3D p2, Point3D p3)
        {
            double[] NormalVector = FaXiangLiang(p1.X, p1.Y, p1.Z, p2.X, p2.Y, p2.Z, p3.X, p3.Y, p3.Z);//先求出法向量
            if (NormalVector[2] < 0)
            {
                NormalVector = new double[3] { -NormalVector[0], -NormalVector[1], -NormalVector[2] };
            }
            Vector3D vector = new Vector3D(NormalVector[0], NormalVector[1],NormalVector[2]);
            double Longitude = 0;
            if (NormalVector[0] > 0&& NormalVector[1] < 0 )
            {
                Longitude = Math.Atan(-vector.X / vector.Y);
            }
            if (NormalVector[0] > 0 && NormalVector[1] > 0)
            {
                Longitude = Math.Atan(vector.Y / vector.X) + Math.PI/2;
            }
            if (NormalVector[0] < 0 && NormalVector[1] > 0)
            {
                Longitude = Math.Atan(-vector.X / vector.Y) + Math.PI;
            }
            if (NormalVector[0] < 0&& NormalVector[1] < 0)
            {
                Longitude = - Math.Atan(vector.X / vector.Y);
            }
            return Longitude * 180 / Math.PI;
        }
    }

    public enum ESatelliteShape
    {
        Sphere,
        SimpleSatellite
    }

    public class Orbit3DParam
    {
        /// <summary>
        /// 轨道半径
        /// </summary>
        public double Radius = 0;
        /// <summary>
        /// 轨道倾角（单位：度）
        /// </summary>
        public double Incline = 0;
        /// <summary>
        /// 升交点赤经（单位：度）
        /// </summary>
        public double Omega = 0;
    }
}
